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Chapter 1 



INTRODUCTION 

1.1 Introduction 

Ordinary differential equations (ODE) and partial differential equations (PDE) have always 
played an important role in modelling physical phenomena. Since only a small number of 
(PDE) are known to have analytical solutions, anyone who wants to develop and use models 
based on (PDE) must be able to obtain a relevant numerical solution efficiently and accurately. 

Numerical methods have played an important role in solving problems with complicated 
geometries and (or) complex governing equations and boundary conditions for which an ana- 
lytical solution is hard or impossible to obtain. Many numerical methods have been developed 
for solving differential and integral equations such as finite-difference methods, finite-element 
methods, boundary element methods, etc. Each of these methods has its inherent advantages 
and disadvantages and the search for easier and more accurate numerical methods is a continu- 
ous and ongoing process. 

The collocation method is a type of numerical techniques used to solve ordinary differen- 
tial equations with boundary or initial conditions. In applying collocation methods, the test 
functions (bases) must be chosen from a complete set to ensure convergence and from lowest 
members of the complete set to ensure high accuracy. An important feature of the colloca- 
tion method, which has contributed to its widespread use, has been the ability to achieve high 
accuracy with few terms in the trial solution. 
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Some of different types of Bases functions used to solve ordinary and partial differential 
equations are Legendre functions, Chebychev functions, Hermite functions, Bessel functions, 
Wavelet basis, Haar wavelet and Sine functions. In this thesis we present the collocation method 
using polynomial bases (B-splines). 

1.2 Layout of dissertation 

The dissertation is organized as follows: 

• Chapter 2 provides a brief outline of the fundamental concepts of B-splines that will make 
the presentation as self-contained as possible. 

• chapter 3, we use the cubic B-spline for solving linear and non-linear second order time 
dependent problem which arises in different physical applications, the error analysis of 
the method is described and this work is published in [1]. 

• chapter 4 we introduce The second order time dependent Emden-Fowler problem which 
has singularity at x — 0, then use cubic B-splines for its solution, this work is published 
in [2]. 

• chapter 5 we use the quintic B-spline for solving linear fourth order time dependent prob- 
lem, the error analysis of the method is described and this work is submitted. 

• Chapter 6 use the quartic B-spline to solve Falkner-Skan equation and we use new algo- 
rithm based on Newton and Secant methods to solve non-linear system. 

• Chapter 7 We give some conclusions and outlines ideas for further research. 

1.3 Work Objective 

The main objective of the work of this thesis can be summarized in the following points . The 
objective has been fulfilled through the completion of the following tasks: 

• Studying the fundamental concepts of the B-spline functions with different orders. 
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• Using B-spline bases with finite difference scheme for solving linear and non-linear sec- 
ond order time-dependent problems. 

• Studying the Emden-Fowler problem and a and present a reliable new algorithm based on 
B-spline to overcome the difficulty of the singular point at x = 0. 

• Using B-spline bases with finite difference scheme for solving linear fourth order time- 
dependent problems. 

• Applying the B-spline methods to obtain solution of Falkner-Skan equation. 

• Deduce complete error analysis for such problems. 
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Chapter 2 
B-SPLINES 



2.1 Introduction 

Spline functions have some attractive properties. Due to the being piecewise polynomial, they 
can be integrated and differentiated easily. Since they have compact support, numerical methods 
in which spline functions are used as a basis function lead to matrix systems including band 
matrices. Such systems have solution algorithms with low computational cost. 

The rapid development of spline functions is primarily due to their great usefulness in appli- 
cations. Classes of spline functions possess many nice structural properties as well as excellent 
approximation powers. 

Splines have many applications in the numerical solution of a variety of problems in applied 
mathematics and engineering. Some of them are, data fitting, function approximation, Integro- 
differential equations, optimal control problems, Computer-Aided Geometric Design (CAGD), 
Wavelets and so on. Programs based on spline functions have found their way in most of 
computer applications. 

2.1.1 Historical note and literature survey 

It is commonly accepted that the first mathematical reference to splines was made in the year 
1946 in an interesting paper by Schoenberg [3] , which is probably the first place that the word 
"spline" is used in connection with smooth and piecewise polynomial approximation. However, 
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the ideas have their roots in the aircraft and shipbuilding industries. Splines are types of curves, 
originally developed for shipbuilding in the days before computer modeling. Naval architects 
needed a way to draw a smooth curve through a set of points. 

The solution was to place metal weights (called knots) at the control points, and bend a 
thin metal or wooden beam (called a spline) through the weights. The physics of the bending 
spline meant that the influence of each weight was greatest at the point of contact and decreased 
smoothly further along the spline. To get more control over a certain region of the spline, the 
draftsman simply added more weights. This scheme had obvious problems with data exchange. 
There was a need for mathematical way to describe the shape of the curve. Cubic Polynomials 
splines are the mathematical equivalent of the draftsman's wooden beam. 

Through the advent of computers, splines have gained more importance. They were first 
used as a replacement for polynomials in interpolation and then as a tool to construct smooth 
and flexible shapes in computer graphics. 

As late as 1960, there were no more than a handful of papers mentioning spline functions 
by name. Some of the papers which have made great contributions in the development of 
splines include (Ahlberg and Nilson [4], Brikhoff and Garabedian [5], Loscalzo and Talbot 
[6], Maclaren [7], Rubin and Khosla [8], Sastry [9], Schoenberg [3]). Univariate splines were 
studied intensely in the 60's, and by the mid-70's they were sufficiently well understood to 
permit a fairly comprehensive treatment in books form. Some of the books which discuss 
splines thoroughly include (Ahlberg et al. [10], de Boor [11], Prenter [12], Schumaker [13], 
Shikin and Plis [14], Spath [15]). 

2.1.2 B-spline for solving differential equations 

The B-spline functions are used as a basis functions in finite element method, Galerkin method 
and collocation method to construct numerical methods for the approximate solutions of B VPs 
occurring in various engineering applications. 

Application of the Galerkin method accompanied by higher-degree polynomials results in 
a higher-degree matrix system. That brings a burden for numerical analysis, and the computa- 
tional cost of the matrix system increases in the evaluation of both linear and nonlinear systems. 
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On the other hand, the collocation method together with B- spline approximations represents 
an economical alternative since it only requires the evaluation of the unknown parameters at the 
grid points. 

Papers on solving ODE's and PDE's using B-splines are quite numerous of which we men- 
tion some with quartic and quintic B-splines. In Refs. (Caglar et al. [16], Saka and Dag [17] 
quartic B-splines are used and solution of various PDE's is obtained by employing quintic B- 
spline collocation method in (Ramadan and El-Danaf [18], Raslan [19], Saka [20], Saka and 
Dag [20], Zaki [21]. 

Some of the earliest papers using spline functions for smooth approximate solution of or- 
dinary and partial differential equations (PDE's) include (Albasiny and Hoskins [22], Bickely 
[23], Crank and Gupta [24], Jain and Aziz [25], Jain and Aziz [26], Rubin and Khosla [8], Sas- 
try [9], Usmani [27], Usmani and Sakai [28], Usmani and Warsi [29]) . These papers demon- 
strate the approximate methods of solving second-, third- and fourth-order linear boundary- 
value problems (BVP's) and solution of elliptic and parabolic equations by spline functions of 
various degrees. 

More recent papers like Quartic spline methods to solve one-dimensional telegraphic equa- 
tions was proposed by Liu [30], second-order one space dimensional hyperbolic equation solved 
by Ding [31], one-dimensional heat equation solved by third degree B-spline functions by 
Caglar 1 [32], nonlinear Burger's equation solved with modified cubic B-splines collocation 
method is proposed by Mittal [33]. 

Today, there are hundreds of research papers on this subject, and it remains an active re- 
search area. 

In the following sections we will introduce some mathematical overview about B-splines 
and how it can be derived from the divided difference table, and some important properties 
which are employed in developing techniques for the solution of differential equations. 
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2.2 The polynomial m splines S m (^) 



Let 7r be a partition t 0 < t± < t n (n > m + 1) of [a, b]. The space S m (ir) of all functions 
from C^" 1-1 ) [a, b] that reduce to polynomials of degree m on each subinterval (ti, t i+ i) of [a, b] 
is called the set of mth-degree polynomial splines with knots at n. The splines of order 1 are 
step functions, those of order 2 are broken lines, and so on. A typical spline of order k is the 
truncated power function with power k. 

Definition 2.1. Basis: A basis for a finite-dimensional space U is a sequence of elements (fi)f =1 
of U such that each / e U has a unique representation / = J2 i= aifi. The number n is the 
dimension of U. 

Theorem 2.2. 77ze monomials and truncated power functions 



form a basis for the spline space S m (ir). In particular, the dimension of this space is k + m — 1. 
proof: See [34]. 

But the truncated power basis has several disadvantages and is not well suited for numerical 
computations. The basis functions are not local; i.e., they are nonzero on the whole interval 
[a, b]. Further, they are almost linearly dependent when the knots are close. Therefore, this 
basis yields dense ill-conditioned linear systems for various tasks. A more suitable basis will 
be introduced in the following section is called B-Splines. 

2.3 Computation of B-spline functions B m (t) 

[12] In this we will introduce a more convenient B-spline basis for the linear space S m (7r) 
of spline functions of degree k. In our work we will be interested in equidistance knots with 
increasing step h where h = and the B-spline functions B m (t) are computed using the m th 
forward difference of truncated power function of degree m as following. 



{1, x, • , (a? •£].)+ ? ' ' ' j (•£ •^m— 1)+ } 




(2.1) 
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where the truncated power function (ti — t)™ is m — 1 continuously differentiable function 
defined as 

{(U -t) m ,wheret < U\ 
(2.2) 
0, wheret > t { . 

But in order to simplify the calculation of B- spline function of higher order Carl de Boor[ll] 
derive the recurrence relation. 

Bi, k (t) = t ~ U + Bi, k ^(x) + U+k ~ X B l+l>k ^(t) (2.3) 

ti+k-l ~ ti ti + k — tj+i 

We start by considering the space Si(n). This consists of piecewise constant functions and a 
basis is. 

{l x e [U,t i+1 ) 

(2.4) 
0 otherwise 

By use of zero degree B-spline function we could easily calculate 1st, 2nd and higher order 
functions. For example if the knot vector is the uniform sequence such that {tj}™ =0 = {i}i=o- 
Here are examples for the first few values of k. For k — 1, the B-spline functions are defined by 

l fe[fi,f i+ i) 
B,At)= { (2.5) 

0 otherwise 

This results in the step functions shown in Figure (2.1) . Notice how each step is closed on the 
left and open on the right and how B a (t) is nonzero only in the interval [t i: t i+1 ) (this interval 
is its support). 

It is also clear that each of them is a shifted version of its predecessor, so we can express 
any of them as a shifted version of the first one and write Bn(t) = B 01 (t — i). For k = 2, the 
functions can be calculated for any % from Equation 

B Q , 2 (t)= t -^-B m (t)+ t -^-B ll (t) 

ti — to £2 — n (2 6) 

= tB 0l (t) + (2-t)B u (t) 

t 0 < t < 1 

\2-t 1 < t < 2 

otherwise 
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) 

Boi(t) 
0 1 



i. J 



2 



0 



B 2l(t) 



ii 



Bai{t) 



Figure 2.1: Uniform B -Spline Functions for k = 1. 
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Bi fl (t) = ^-TB 11 (t) + p-^-B 21 (t) 

l 2 — z l z 3 ~~ l 2 

= (t-l)B n (t) + (3-t)B 21 (t) 



#2,2 (*) 



t-1, 

3-f, 
0, 

t-t 2 



when 1 < t < 2 
when 2 < t < 3 
when otherwise 

^3 — r 2 C4 — ^3 

= (t - 2)B 21 (t) + (4 - t)B 31 (t) 

t - 2, when 2 < t < 3 

4 - t, when 3 < t < 4 

0, when otherwise 



(2.7) 



(2.8) 



The hat-shaped functions are shown in Figure (2.2). Notice how B i2 (t) spans the interval 
\ti,t i+2 ). It is also obvious that each of them is a shifted version of its predecessor, so we 
can express any of them as a shifted version of the first one and write B i2 (t) = B 02 (t — i). 



For k = 3, the calculations are similar 



£2 — H) £3 — H 

^ ^ / \ 3 t „ , . 
2^o 2 (t) + — B 12 (f) 



0 < t < 1 



5(2-t) + ^(2-l) 



1< i < 2 



(3-t) 2 
2 

0 

5l )3 (t) 



2 < i < 3 
otherwise 

^B 12 (t) + ^B 22 (t) 
-s-Si 2 (i) + -— fl 22 (t) 



(2.9) 



(2.10) 
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Figure 2.3: Uniform B-Spline Functions for k = 3 



1 < t < 2 



2 

2 < t < 3 



-2f 2 + 10f-ll 



2 

3 < t < 4 



2 

0 otherwise 



In general, the support of B ik {t) is the interval [tj, and B ik (t) = B 0k (t — i). Figure (2.4) 
shows how a general function B ik {t) is constructed recursively. Each B ik {t) function in this 
triangle is constructed as a weighted sum of the two functions immediately to its left. 
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B j.fc-2 



B 



B 



B 



i.k 



B 



B 



B i+Jfe-2,1 



Figure 2.4: Triangular array for computing the B-splines of order k — 1, 2, 3, ■ • • 



2.4 Why B-splines 



B-splines are mathematically more sophisticated than other types of splines (M-spline, L-spline, 
etc.)- The B-spline curve has superior properties that make them suitable for shape representa- 
tion and analysis. Some of the important properties are: [35] 

2.4.1 Smoothness and continuity: 

Due the differentiability properties of the basis B-spline function, the B-spline curve is at least 
continuously differentiable of up to p — 1 times. For the cubic B-spline curve, for example, the 
2nd derivative B-spline curve exists. 

2.4.2 Boundedness: 

The B-spline curve is confined in the convex hull of its control polygons. In fact, for t in the 
knot span [tj, t i+1 ), the B(t) is bounded by the convex hull of the control points Pi- P , ■ ■ ■ ,Pi. 
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2.4.3 Local shape controllability: 

Due to the local support of the basis B-spline function, moving p, L affects r(t) only in the interval 
[U, U+p+i) and hence the change of shape produced by varying one of pt is locally confined. 

2.5 Properties of B -splines 

The B-splines of order k have the following properties: 

2.5.1 Spline structure. 

B-spline Bi is indeed a spline, i.e., a piecewise polynomial function of order k with the knots 

ti, tj+fc.-Bj G C k 1 m 

2.5.2 Finite support. 

B-splines have a finite support 

Supp^ = [ti,t i+k ] 

2.5.3 Peano kernel for divided difference. 

B-spline Bj is the Peano kernel in the integral representation of the divided difference functional 

[ti...ti + k\- 

2.5.4 Positivity. 

Since B^\ is just a step function with support on [t i: t i+ i], namely 

{l x e [U,t i+1 ) 

0 otherwise 

and the recurrence (2.3) provides immediately finiteness of support of B-splines and their posi- 
tivity 

Suppfii = [ti,ti+k] B itk >0 on (ti,ti+k) 
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SOLUTION OF SECOND ORDER TIME 



DEPENDENT PROBLEM 



3.1 Introduction 

In this chapter, we employ the B-spline method to solve linear and non-linear 2 nd order time 
dependent partial differential equations in the form 

— + {a- l)A(x,t) — = u xx + (a - 2)B(x,t)u x + C(x,t)u + f(x,t,u,u t ) (3.1) 
subject to boundary conditions 

u(p,t)= gi (t), u(l,t)=g 2 (t) (3.2) 

and initial conditions 

u(x, 0) = fJ>i(x) 

(3.3) 

(a - l)u t (x, 0) = (a - 1)/jl 2 (x) 

where a G {1,2} and pLi,^ are prescribed functions and A, B and C, are finite constants 
or continuous functions. The case a = 1 corresponding to diffusion-convection problem for 
C — 0, reaction problem for B = 0 and heat problem for B = C = 0, respectively. The case 
a = 2 corresponds to telegraph problem and wave problem for B = 0. respectively. Several 
numerical methods have been developed for solving special types of time-dependent problem 
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such as finite-difference [36], finite volume [37], finite-element [38], boundary-element meth- 
ods [39], sinc-galerkin [40], wavelet-galerkin [41, 42], multi-grid methods [43], the decomposi- 
tion method [44] , the dual reciprocity boundary integral equation [45], the Chebyshev cardinal 
functions [46], spline methods [47, 48, 49, 50, 51], etc. The organization of this chapter is as 
follows. Section (3.2), contains notation, definitions and some results of the B-spline. The so- 
lution method is introduced in section (3.3) and (3.4) where four main theorems are presented . 
In Section (3.3), the error analysis of the method is described in detail. Section (3.4) is devoted 
to derivation of the discrete system. Finally, some numerical examples are presented in Section 
(3.5), where the exact solution is explicitly given. 



3.2 Description of B-spline method 



The third-degree B-spline is used to construct numerical solutions to a given problem (3.1). A 
detailed description of third- degree B-spline functions can be found in [52]. The third-degree 
B- splines are defined as : 



B i . 1 (x) = B 0 (x-(i-l)h), 



* = 2,3,.... 



and 



B 0 (x) 



x 



6h 3 



-3x 3 + 12hx 2 - 12h 2 x + Ah 3 , 
3x 3 - 2Ahx 2 + 60h 2 x - AAh 3 , 
-x 3 + Ylhx 2 - A8h 2 x + 6Ah 3 , 



0 < x < h 
h < x < 2h 
2h<x <3h 
3h < x < Ah 



To solve second-order boundary value problems B, B', B" evaluated at the nodal points are 
needed. 

i j 



B i \Xj / 



4, if 
1, if 
0, if 



i-j = ±l 
i-j = ±2 
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Similarly, we can show that: 

0, if; i = j 

if; i-j = ±l , B i( x i) = < 

if; i-j = ±2 



±2 



.12 

/j2 , I J , 



o, 



_6_ 
ft 2 ' 



o, 



I = J 

if; i-j = ±l 
if; i- J = ±2 



The cubic B-spline interpolation is a linear combination of the cubic B-spline basis as follows: 



n-l 



S(x)= Y^CiBiix) 



i=-3 



where Cj are unknown real coefficients and -Bj(x) are third-order B-spline functions 

3.3 Error analysis 
3.3.1 Governing equation 

Let f(x, t, u, u t ) = ij){x, t) we consider the linear non-homogeneous second-order time-dependent 
equation: 



Setting 



0 a u du 

— + (a - l)A(x, t)-^ = u xx + (a - 2)B{x, t)u x + C(x, t)u + ^(x, t) (3.4) 



d a u (a — 1) u i+ i — aui + 



dt a (- A t) a 

Equation (3.4) may be approximated by the following scheme 



(a — l)u i+1 — aui + Ui-i 
(-At) a 



+(a-l)A(x,t i+1 ) 



Ui - Ui-i 



At 



d 2 Ui 
dx 2 



+ (a-2)B(x,t i+1 ) 



du i+ i 
dx 



C(x,t i+1 ) 



--u't +1 +(a-2)B(x,t i+1 )u' i+1 
+ C(x, t i+1 )u i+1 + ip(x, t i+1 ) 

(a-iy 



-Aty 



u i+ i = M(x,t i+1 ) (3.5) 



where 



M(x,t i+ i) = -if)(x,t i+ i) + 



(a — 1) , . . a 
'■A(x,t i+1 ) - 



At 



Ui + 



-Aty 

1 («-l) 
_(-At) a A* 



A(x,t i+1 ) 



Ui-!. 
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Suppose S(x) is the cubic B-spline interpolating u{xk, U+i) at the (i + l) th time level then 



n-l 



S(x) = CkB k (x) 



(3.6) 



k=-3 



Lemma 3.1. Let the collocation approximation S(x) to the solution u(x) of the boundary value 
problem (3.5) be (3.6), then the following relation holds 



S ^xy^i) ^ gg Uxxxxx (pCi) 0(H ), 

h 2 h 4 
12 36U 



(xi) + 0(h 6 ) 



Proof. Equation (3.6) can be simplified to 



u(xi) S(xi) 

and we can easily get 

u x( x i) ~ S (Xi) = ( 
Uxx{%i) ~ S {Xj^j 

the following relations can be obtained [48] 

/< ^'(^-i) + \s\ Xi ) + l -S\x i+1 ) 

6 3 6 



Ci-l + 



l) Ci+ (l ] Q+1 



1 



Ci-l 



1 

2/1 

2 



Ci+l 



W) Ci+ {h? ]Ci+1 



= h 




IK 







-Q_ 2 + -Q ] 4 



2 /-l 1 \ 1 f-1 1 

3 l x 2^ Q - 1 + 2h Cl+1 ) + 6 V 2ft 04 + 2h Cl+2 



-1 2 2 1 

6 3 3 6 



Adding and subtracting \ci we get 

h 



^S"(a;i-i) + ^'(xi) + ]:S'(x i+1 ) 
6 3 6 



-1 2 112 1 

6 3 6 6 3 6 



which reduced to the following relation 

h 



(3.7) 
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Similarly we can find formula for S"(xi) as following 



h 2 S"( Xi ) = h 2 



h 2 S"(x t ) = 6 



" 1 2 1 

y 2 Ci-i - + - 

Ci + 4q+i + c i+2 \ 



... , -i-i + 4cf + Cf+i \ 

e ; v e ; 

(ih + 2^+ 2 ) + 2 (i/^ 1 + 27^ +1 



(3.8) 



2h 



we get the following relation 



h 2 S"( Xi ) = 6 [S(x i+1 ) - S( Xi )] - 2h [S\x i+1 ) + 2S'(x i )\ 
By using the operator notation E(S(xi)) = S(x i+ i) equation(3.7) may be written as [9] 



(3.9) 



h 



1 , 2 1 „ 
+ - + -£ 
6 3 6 



S'( Xi ) = \(E- E- 1 ) u{ Xi ) 



since E = e hD where D = -j- substitute in the above relation 

ax 



h 



l + l(e- hD + e hD ) 



S'(x t ) = - [e 



hD _ e -hD 



] S( Xi ) 



Expand in Maclaurine series in power of hD 

1 fh 2 D 2 h A D 4 h 6 D 6 



h 



1 + 



3 

S\ Xi ) = 
S'(xi) ■ 



2! 



4! 



6! 



S'( Xi ) 



' „ h 3 D 3 h 5 D 5 
hD + —— + 



h 7 D 7 
~7\~ 



3! 



5! 



, h 2 D 2 h A D A h 6 D 6 



6 



72 2160 



1 T h 2 D 3 h A D 5 h 6 D 7 
D + —^ + —r- + —^ + 



3! 



5! 



7! 



h 2 D 3 h A D 5 h 6 D 7 
D + —— + —^ + —^ + 



3! 



5! 



7! 



1 - 



h 2 D 2 h A D A h e D e 



6 



72 2160 



u 



u[x. 



+ 



h 4 D A h 8 D 8 h 6 D G 
36 + 5184 + 216 + 



+ 



where Du(x) = u x (x) and D 5 u(x) = u xxxxx (x) and so on 

•S* (^Cj) Uxi^Xi) Y^^u xxxxx (xi) + 0(/l ) 

Similarly we substitute in (3.9) 

h 2 S"(xi) = 6 [E - 1] Sfo) - 2/i [2 + E] S'{xi) 

= 6 [e hD - 1] S(^) - 2h [2 + e^] S'(^) 



u(xi) (3.10) 



27 



CHAPTER 3. SOLUTION OF SECOND ORDER TIME DEPENDENT PROBLEM 



Expand in power of hD and substitute by S'(xi) 



h 2 S"( Xi ) 



h 2 D 2 h 3 D 



6 -1 + 1 + HD + — — + —— + ... )-2h\ 2 + l + hD + —— + 



2! 3! 



h 2 D 2 h 3 D z 



2! 3! 



„ h A D 5 h 6 D 7 
D-—— + —— + 



180 1512 



u{xi) (3.11) 



Hence we can get 



S"( Xi ) = 



D 



h 2 D A h 4 D 6 



+ + ••• 



12 360 



u(xi) 



Finally 



S"(xi) 



h 2 h A 6 

^xxy^i) ^i^xxxxy^i) T^^^xxxxxxy^i) 0{h ) 



(3.12) 

□ 



Theorem 3.1. The truncation error of time dependent equation is given by 

((« " 2) (At) u t + a { -^u tt + (a- 2) { -^fu m + O ((At) 4 )) 
+A(a - 1) (u t - { -^-u tt + { -^fu ttt + O ((At) 3 )) = X (x,t), 



where 



x{x, i) = u xx + B(a — 2)u x + Cu + ^(x, t) 



Proof: 



(3.13) 



— + A(a- 1)— =w ;c;c + J B(a-2)M :r + C , M + ^(x,t) 



by applying difference formula 

(a — l)Mj+i — cnUi + 
_ (-At)° 

Applying Taylor series expansion 

1 



+ A(a - 1) 



At 



(-At)< 



(At) 2 (At) 3 
a-l)[u+ (At)u t + -^-u tt + ^-u ttt + 



(At) 2 (At) 3 
u - (At)u t + ^-u tt - ^-um + 



.{a — 1) . . . . 

+ ^4 - — (u — it + (At)u t 



At 



{At) 2 {Atf 
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(3.14) 



Rearranging the equation leads to : 

((« " 2)(Af) u t + a { -^u tt + (a- 2) { -^fu ttt + 0 ((At) 4 )) 

+A(a - 1) (u t - { -^-u tt + ( -^u ttt + 0 ((At) 3 )) = X (x,t) 

Theorem 3.2. By using the combination of the finite difference approximation and cubic B 
spline method, the truncation error for the time-dependent equation (3.1) is 

O ((At) + h 2 ) □ 

Proof: Combine Lemma (3.1) and Theorem (3.1). For Eq. (5.24) 

((« - 2) (At) u t + a^u a + (a- 2)^u ttt + O ((At) 4 )) 

+A(a - 1) (u t - { -^-u tt + { -^fu m + O ((At) 3 )) = X (x,t) 
we have two cases: 

Case (i) a = 1 

X(x, t) = u XI - Bm^ + Cu + ^(z, t) 

and equation (3.14) will be 

-1 / (At) 2 (At) 3 \ 

— I (-At) «t + -^-u tt - -^-u ttt + O ((At) 4 ) J = u xx - Bu x + Cu + tP(x,t) 
from lemma (3.1) 

u t - ^u tt + { -^fu ttt + O ((At) 4 ) 



h 2 /i 4 6 

U xx (Xl) ~7~^U XX xx(%i) „ „ „ U XXXXX x 0(h ) 

12 doU 



/i 4 

-B(u x (xi) + — w OT:c (xi) + 0(/i )) + Cm + ^(x, t) 



the truncation error can be estimated as e = u t — x(x, t) 

e = -y-uu — Uut+O ((At) )-— u xxxx (x i ) + —u xxxxxx {x i )+B — 

so the truncation error is 0(At, h 2 ) 
Case (ii) a = 2 



X(x, t) = u xx + Cu + tp(x, t) 
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and equation (3.14) will be 

{{M f Ua + o ((At) 4 )) + A (u t - { -^u tt + ^ ttt + 0 ((At) 3 )) = 

u xx + Cu + ip(x,t) 

which is reduced to 

u tt + 0 ((At) 4 ) + A«t - A { -^u tt + A { -^fu ta + O ((At) 3 ) = 
/i 2 /i 4 6 

the truncation error can be estimated as e = (u tt + Au t ) — t) 
so the truncation error is 0(At, /i 2 ). 



3.4 Scheme of the method 
3.4.1 Linear time-dependent 



Let the solution u(x, t i+1 ) of the problem (3.5) be approximated by 



n-1 

u(x,t i+1 ) = ^ CjBj(x) (3.15) 

j=-3 



where c,- are unknown real coefficients and Bj(x) are the third-degree B-spline functions. Let 

x 0 , xi, x n be n + 1 grid points in the interval [a, b] so that = a + ih,i = 0, 1, 2, n 

where 



x 0 = a, x n = b and /i = 

n 



and 

n— 1 rt— 1 

J=-3 i=-3 

Theorem 3.3. If the assumed approximate solution of the problem (3.5) is (3.15) then the dis- 
crete collocation system for the determination of the unknown coefficients {Cj}"~^ 3 is given 
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by 



n—1 n—1 

J2 c j B'j(x k ) + (a-2)B(x k ,t i+1 ) Cj B j(^k) + 

J=-3 j=-3 



C(x k ,t i+1 ) 



(a -1) 



n-1 



C 3 B Mk) 

J=-3 



-At) a _ 

M(x k ,t i+1 ), fc = 0,l,...,ra (3.17) 



and boundary conditions (3.2) can be written as 



n-1 



i=-3 

n-1 



^ CjBj(x) = g2(t), for x = 1. 

j=-3 



□ 



Proof: If we replace each term of (3.5) with its corresponding approximation given by 
(3.15) and (3.16) and substituting x = x k and applying the collocation to it.D 
Hence, the initial condition for i = 0 we have 

uq = u(x, 0) 
Then the system in (3.17) takes the matrix form 



A C = I 



(3.18) 



where 



C = 



/ \ 

C-3 

C-2 



C n -2 

y c n -i J 



( 



F = 



9i(t) 
M(x 0 ,t) 
M(x u t) 



M(x n _!,t) 
M(x n ,t) 
V 92{t) J 
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also the matrix A can be written as 



A = 



14 10 

Si(x 0 ) S 2 (x 0 ) 5 3 (x 0 ) 0 

0 <5l(Zi) S 2 (Xi) S 3 (x 1 ) 



V 



0 
0 



0 
0 



0 
0 



0 
0 



0 0 

0 0 

0 0 

8i(x n ) 5 2 (x n ) S 3 (x n ) 

1 4 1 



Also the coefficients of the matrix A have the following form 



Si{x k ) 
S 2 (x k ) 
5 3 (x k ) 



^)+(a-2)B(x k ,t i+1 ) ) + 



C(xk,t i+ i) 



-12 

6 



+ 



C{Xk, t 



(-Af)° 
■3 



•4, 



} . 2 ; + (a -2)B(x k ,t i+1 ) , 



- C(x k ,t i+1 ) 

Theorem 3.4. The matrix A in the system (3.18) is nonsingular. 



(q-i) 

(-At)" 



(<*-!) ' 
(-At)" 



Proo/ The basis function B^x) is defined only in the interval and outside of 

this interval it is zero. Therefore,^ (x) is having non- vanishing values at the mesh points 
{xj_i,Xj,x i+ i}and at other mesh points the value of Bi(x) is zero. It is clear that from the 
definition of cubic B-splines defined in section 3.2, the value of B^x) at Xj is dominating for 
j = i when compared with the values of Bi(x) when j ^ i. The derivatives of B^x) up 
to second order also have the same nature at the mesh points as in the case of Bj (x) . Using 
these facts, we can say that the matrix A is a tri-diagonal band matrix with nonzero entries and 
dominant principal diagonal elements. Hence the matrix is nonsingular. □ 



Now we have a linear system of n + 3 equations for the n + 3 unknown coefficients, namely, 



n— 1 



We can obtain the coefficients of the approximate solution by solving this linear 



system by Q-R method. 



3.4.2 Non-linear time-dependent 
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We consider the non-linear non-homogeneous time dependent equation (3.1)and by apply- 
ing finite difference, 



d 2 Ui+i , , n \ -n, . xdu i+ i 



dx 2 
where 



+ (a-2)B(x,ti 



) H+lj 



dx 



+ 



C(x,t i+1 ) - 



(Q-l) 
{-At) a 



u. 



+1 



N(x,t 



i+lj 



(3.19) 



iV(a;,f i+ i) = -/ ( x,t i+1 ,Ui, — ^ 1 | 



(a — 1) a 
-A(x,t i+1 ) - 



At 



1 



(-At) a 
a -I) 



A(x,t 



i+l) 



Ui-1- 



( — At) a At 

Theorem 3.5. If the assumed approximate solution of the problem (3.19) is (3.15) then the 
discrete collocation system for the determination of the unknown coefficients {Cj}™I^ 3 is given 



by 



n—1 n—1 

J2 CjB'^Xk) + (a - 2) B(x k ,t t+l ) c j B j( x k)+ 

J=-3 j=-3 

n-l 

Y c jBj(x k ) = N(x k , t i+l ), k = 0, 1, n. (3.20) 

J=-3 



C(x k , t i+ i) — 



(a -I) 



( — At) a 

Proof: If we replace each term of (3.19) with its corresponding approximation given by 
(3.15)and (3.16) and substituting x = x k and applying the collocation to it.D 

Now we have a non-linear system of n + 3 equations in the n + 3 unknown coefficients. 
We can obtain the coefficients in the approximate solution by solving this non-linear system by 
Newton 's method [53] . 
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3.4.3 Newton's Method 



To solve the system of equations (3.20), We write it in the form 



Fi(c_ 3 ,C_2, ...,c„_i) 0 



F(C) = 



F 2 (C_ 3 ,C_2, ...,C n _i) 



-Pn+3( c -3 ? c -2, C n _i) 



V 



/ V / 



(3.21) 



where C is the column vector of independent variables and F is the column vector of the func- 
tions Fj, with Fj(C) = Fj(c_3, c_ 2 , c n _i), 1 < i < n + 3. The number of functions that are 
set to zeros is equal to the number of independent variables. A very good method for solving 
equation (3.21) is Newton method. 

Let be the guess at the solution for iteration j. Let F^ denote the value of F at the j th 
iteration. Assuming that || F^ || is not too small, we seek update vectors Ac^ 



,C?+i) = c 0") + Ac (i) 



c j+1 



J+i 

c -2 



C> 



Ac> 



Ac 7 



\ 



(3.22) 



V 4-1 ) \ 4-i / V A 4-i / 

such that F(c J+1 ) = 0. Using the multidimensional extension of Taylor's theorem to ap- 
proximate the variation of F(c) in the neighborhood of c l gives : 



F(c J + Ad) = F(c (i) ) + F'(c (j) )Ac (j) + o(|| Ac (j) f) 



(3.23) 
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where F'(c^) is the jacobian of the system of equations. 



/ dFi 



F'(c)=J(c) = 



dc- 



dF 2 
dc-i 



dc- 



W lft(c) 



dF n+3 f„\ 8F n+3 



He) ^(c 



V 



dFi 

dc n - 



dF 2 
dc n - 



\ 



/ 



(3.24) 



Neglecting higher order terms and designating as the Jacobian evaluated at c^, we can 
rearrange equation (3.23). 



F(c w + Ac (i) ) = F(c (j) ) + J (j) Ac 



(3.25) 



The goal of Newton iterations is to make F(c^ + Ac^) = 0 , so setting that term to zero in the 
preceding equation gives 



j(i) Ac (i) = _ F ( c (i)) (3.26) 

Equation (3.26) is system of n + 3 linear equations in the n + 3 unknown Ac^\ Each Newton 
iteration step involves evaluation of the vector the matrix and the solution to equation 
(3.26). A common numerical practice is to stop the Newton iteration whenever the distance 
between two iterates is less than a given tolerance, i.e., when || c^ +1 ^ — ||< e. 



3.4.4 The initial state 

The initial vector can be determined from the initial condition u(x, 0) = pLi(x) which gives 
n + 1 equation in n + 3 unknowns [51]. For the determination of the unknowns we need two 
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more conditions can be determined as follows: 



« x (0,0) =/x' 1 (0) 



^(1,0) =^(1) 



the initial vector is then determined as the solution of the matrix equation 
/ 



" 0 
4 



0 1 



3 
h 

1 

4 



0 
0 
1 



0 
0 



0 
0 
0 



C-2 



0 0 0 



1 4 1 



V 



o o 0 ••• f o I j 



Cn-2 

y c n -i j 



A*i(0) 



/ii(l) 



(3.27) 



3.5 Examples and comparisons 



We have implemented our method on five examples, three linear and two nonlinear which were 
collected in [40, 54, 55, 56, 57]. The performance of the B-spline method is measured by the 
maximum absolute error e k which is defined as : 

&k Inexact UBspline] 

Example 3.1. [54] In the case a = 1, B = 0,C = — (1 + Ax 2 ) equation (3.1) becomes 



ut = u xx ~ (1 + Ax 2 )u 



subject to boundary conditions 



l+t 



u(0,t) = e*, u(l,t) = e 
with the initial condition 

u(x,0) = e x2 

whose exact solution is 

u(x,t) = e t+x ' 2 

The Maximum absolute error is displayed in Table (3.1) for different n at different times with 
At = 0.001 . The numerical results for B-spline solution, exact solution and absolute error are 
illustrated in Figure 3.1 
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n 


T = 0.1 


T = 0.5 


T= 1 


10 


3.05 E-03 


5.79 E-03 


9.56 E-03 


20 


7.18 E-04 


1.37 E-03 


2.26 E-03 


40 


1.44 E-04 


2.68 E-04 


4.43 E-04 


60 


3.91 E-05 


6.75 E-04 


1.11 E-04 


75 


1.14 E-05 


1.89 E-05 


3.13 E-05 



Table 3.1: Results for example 3.1. 
Example 3.2. [56] In the case a = 2, A = 20, C = —ft 2 equation (3.1) becomes 

u tt + 20u t + f3 2 u = u xx + [3 + /3 2 - 40]e~ 2t sinh(:r) 
subject to the boundary conditions 

«(0,t) = 0, u(l,t) = ^(e-e' 1 )e- 2t , 

with initial conditions 

u(x,0) = sinhx, u t = — 2sinhx, 

whose exact solution is 

u(x,t) = e~ 2 *sinhx. 

In Table (3.2), we compare the results obtained by the B-spline method at time t = 1.0 with 
those obtained by Mohanty [56]. The numerical results for B-spline solution, exact solution and 
absolute error are illustrated in Figure 3.2 



n= 16 


B-spline 


Results by [56] 


0 = 10, £ = 5 


6.44 E-04 


6.752 E-04 


e = 20, p = 10 


4.44 E-04 


4.496 E-03 


0 = 40,13 = 4 


1.714 E-03 


1.866 E-03 



Table 3.2: Results for example 3.2. 
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approximate solution 




Exact solution 






(a) Approx. solution 



(b) Exact solution 



Absolute error 



x 10 



1.5- 



0.5- 




(c) Absolute error 



Figure 3.1: Results for example 3.1 at At = 0.01 and n = 75 



Example 3.3. [40] In the case a — 2, A — 1, C — — 1 equation (3.1) becomes 

Utt + u t + u = u xx + [2 cos 2t + sin 2t + sin 2 1] (x 2 — x 3 ) + (6x — 2) sin 2 1 



subject to the boundary conditions 



with initial conditions 



u(0,t) = 0, «(!,£) = 0 



u(x,0) = 0, u t (x,0) = 0 
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approximate solution 



Exact solution 



15 i 




0.5- 



1.5 i 



0.5- 



0 0 




(a) Approx. solution 



(b) Exact solution 



Absolute error 




(c) Absolute error 

Figure 3.2: Results for example 3.2 at At = 0.01 , n = 16 , 9 = 10 and 0 = 5 
whose exact solution is 

u(x, t) = x 2 (l — x) sin 2 1 

The Maximum absolute error is displayed in Table (3.3) at times T = 0.01, T = 0.5 and 
T = 1.0 with At = 0.001 . The numerical results for B-spline solution, exact solution and 
absolute error are illustrated in Figure 3.3 

Example 3.4. [55] In the case a = 2,A = C = 0 and f(x,t,u,u t ) = ^(u 2 — l)u t + [tt 2 + 
7 2 e -27 * sin 2 7ix]e~" /t sin(7rx) equation (3.1)reduced to the Van der Pol type non-linear wave 
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X 


T = 0.01 


T = 0.5 


T= 1.0 


0.1 


8.879 E-08 


5.799 E-06 


4.885 E-06 


0.2 


3.144 E-07 


1.104 E-05 


1.132 E-06 


0.3 


6.188 E-07 


1.518 E-05 


1.936 E-06 


0.4 


9.428 E-07 


1.738 E-05 


2.759 E-06 


0.5 


1.228 E-06 


1.706 E-05 


3.427 E-06 


0.6 


1.414 E-06 


1.487 E-05 


3.769 E-06 


0.7 


1.443 E-06 


1.181 E-05 


3.631 E-06 


0.8 


1.257 E-06 


8.262 E-06 


2.918 E-06 


0.9 


7.946 E-07 


4.277 E-06 


1.644 E-06 



Table 3.3: Results for example 3.3. 



equation 

Utt = Uxx + l(u 2 — l)u t + [tt 2 + 7 2 e~ 27 * sin 2 nx]e~ jt sin(7nr) 
subject to the boundary conditions 

u(0,i) = 0, u(l,f) = 0 

with initial conditions 

u(x,0) = sin7nr, u t (x,0) = — 7sin7rx 

whose exact solution is 

u(x,t) = sin(7ra;)e _7 *. 

The maximum absolute error is displayed in Table (3.4) for 7 = 1 at times T = 0.01, T = 0.1 
and T = 0.2 . 

Example 3.5. [57] In the case a = 2,A = C = 0 equation (3.1) becomes 

u u = ~u xx + 4 u 3 

subject to the boundary conditions 

u(0,f) = — , u(l,t) = 
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approximate solution for Ex(1) 



Exact solution for Ex(1) 



0.15 




0 0 




(a) Approx. solution 



(b) Exact solution 



Absolute error for Ex(1) 




(c) Absolute error 



Figure 3.3: Results for example 3.3 at At = 0.01 and n = 10 



with initial conditions 



u(x, 0) 



whose exact solution is 



1 

4 + x' 



u(x, t) 



u t {x,0) 



(4 + x) 



1 



4 + x + t 

The maximum absolute error is displayed in Table(3.5) at times T = 0.01, T = 0.1 and T = 0.2 
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X 


T 1 A A 1 

1 = 0.01 


T A 1 

1 = 0.1 


T 1 A O 

1 = 0.2 


0.1 


r\ s s A T! f\s- 

2.661 E-06 


1 O A A T" J A A 

1.349 E-04 


A AAAC 1 

0.0005 1 


0.2 


5.091 E-06 


2.771 E-04 


A AA 1 OA 

0.00120 


0.3 


7.059 E-06 


A 1 /CO T7 A A 

4.162 E-04 


A A A 1 /C A 

0.00169 


0.4 


8.347 E-06 


5.224 E-04 


A AAO 1 A 

0.00219 


U.j 


0. /yt) li-UD 


C S>^C TJ f\/i 

J.OZJ I1-U4 


U.UUzjo 


0.6 


8.347 E-06 


5.224 E-04 


0.00219 


0.7 


7.059 E-06 


4.162 E-04 


0.00169 


0.8 


5.091 E-06 


2.771 E-04 


0.00109 


0.9 


2.661 E-06 


1.349 E-04 


0.00052 



Table 3.4: Results for example 3.4. 



X 


T = 0.01 


T = 0.1 


T = 0.2 


0.1 


1.234 E-07 


2.898 E-06 


2.383 E-05 


0.2 


1.291 E-07 


5.687 E-07 


1.843 E-05 


0.3 


1.166 E-07 


1.657 E-06 


1.805 E-05 


0.4 


1.102 E-07 


1.019 E-06 


7.662 E-06 


0.5 


1.024 E-07 


1.214 E-06 


9.448 E-06 


0.6 


9.812 E-08 


9.379 E-07 


4.781 E-06 


0.7 


8.455 E-08 


1.183 E-06 


1.074 E-05 


0.8 


1.071 E-07 


5.072 E-07 


9.378 E-06 


0.9 


1.234 E-07 


1.601 E-06 


1.262 E-05 



Table 3.5: Results for example 3.5. 
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Chapter 4 



NUMERICAL SOLUTIONS FOR THE 



TIME-DEPENDENT EMDEN-FOWLER 



TYPE 



4.1 Time-dependent Emden-Fowler Equation 

In this chapter, we employ the B-spline method to solve the time dependent partial differential 
equations. For clarity, the method is presented for the heat equation 

T 

Ut = u xx H — u x + ai](x,t,u) + fi(x,t), 0 < x < 1, t>0 (4.1) 

Ob 

subject to boundary conditions 

u x (0,t) = 0, *)=£(*), (4.2) 

and the initial condition 

u(x,0) = u(x). (4.3) 

Some forms of the above equations model several phenomena in mathematical physics and 
astrophysics such as the diffusion of heat perpendicular to the surface of parallel planes, theory 
of stellar structure, the thermal behavior of a spherical cloud of gas, isothermal gas sphere and 
theory of thermionic currents [58, 59, 60]. 
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The solution of the time dependent Emden-Fowler equation as well as variety of linear 
and nonlinear singular IVPs in quantum mechanics and astrophysics is numerically challenging 
because of singularity behavior at the origin. 

The singularity behavior that occurs at the point x = 0 is the main difficulty in the analysis 
of equation (4.1)-(4.3). The approximate analytical solutions to several forms of the above 
problems were presented by Shawagfeh [61] and Wazwaz [62, 63, 64, 65] using the Adomian 
decomposition method (ADM) [11]. Very recently, Chowdhury and Hashim [66] and Belal 
et al. [12] respectively applied the homotopy -perturbation method (HPM) and the variational 
iteration method (VIM) [13] to obtain approximate analytical solutions of the time-dependent 
EmdenFowler type equations. 

The chapter is organized as follows. Section (4.2) we introduce error analysis for proposed 
problem, we formulate our B-spline collocation method to the solution of (4.1)-(4.3) in section 
(4.3). In Section (4.4), numerical experiments are tested to demonstrate the viability of the 
proposed method. 



4.2 Error analysis 
4.2.1 Governing equation 

Since x = 0 is singular point of equation (4.1), we first modify equation (4.1) at x = 0. By 
L'Hospital rule [32], the boundary value problem (4.1) is transform into 

u t = p(x) u xx + q(x) u x + ar)(x, t, u) + p,(x, t), 0 < x < 1, t>0 (4.4) 



where 

r + 1, x=0; 0, x=0; 

p(x) = { q(x) = < 

1, x^O ( r -, x^O 

Let 7](x, t, u) = h(x, t) u, we consider the linear non-homogeneous Emden-Fowler differential 
equation, difference schemes for this problem considered as following [16, 50]: 



p(x] 



d 2 u i+x du i+x 

, 2 +OT — j — 

ax z ax 



+ a h(x, t) u i+ i + fi(x, t) = 



At 



(4.5) 
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or 



(A t) p(x) u\% + A t q(x) u\% + [a (A t) h(x, t) - 1] u i+l = - Ui - (A t) fi{x, t) 
Suppose S(x) is the cubic B-spline interpolating u(x k) t i+1 ) at the (i + l) th time level then 

n-l 

S(x) = C k B k (x) (4.6) 

fc=-3 

Theorem 4.1. [48] Let the collocation approximation S(x) to the solution u(x) of the boundary 
value problem (4.1) is (4.6), then the following relation is hold. 

h 4 



h 2 , , h 4 



loll 



S [%i) U xx (Xi) -^Uxxxxi^i) ~\~ ggg ^xxxxxx 0(ll ) 

4.2.2 Truncation error of time dependent Emden-Fowler equation 



Theorem 4.2. By using the combination of the finite difference approximation and cubic B- 
spline method, the truncation error is O (At, h 2 ) 

proof : Consider the equation (4.5) when apply finite difference 
1 

— (u i+ i - = p(x)u xx + q(x)u x + ah(x, t)u + /i(x, t) 

apply Taylor expansion for L.H.S 

1 , * (At) 2 (At) 3 . . , , . . 

— (Mi + Atu t H — u tt H ^— u ttt H = p{x)uxx + q{x)u x + ah{x, t)u + «(x, t) 

u t + ~^r u tt + 0 ) u ttt + • • • = p{x) ( u xx — -rzh 2 u XX xx + 7^.h 4 u xxxxxx + 0(h e 



2! 3! V 12 360 

+q(x) (u x - ^h A u xxxx + 0(/i 6 )^ + ah(x, t)u + «(x, t) 

if 

p{x)u xx + q(x)u x + ah(x, t)u + p,(x, t) = X(x, t) 
we can calculate the truncation error which is defined as = u t — \(x, t) as: 
At (At) 2 (h 2 h 4 \ rt/Ifi , 

6 = ~¥. Utt + ~^T Uttt + " ' + V 12 ~ 180 ) Uxxxx ~ ' ^ 



46 



CHAPTER 4. NUMERICAL SOLUTIONS FOR THE TIME-DEPENDENT EMDEN -FOWLER TYPE 



4.3 B -spline solutions for linear time-dependent Emden-Fowler 
4.3.1 Linear time-dependent 

Let 

n-l 

«(*)= ^CjBiix) (4.7) 

J=-3 

be an approximate solution of equation (4.1), where Q are unknown real coefficients and 

are third-degree B-spline functions. Let x 0: x 1: . . . , x n be n + 1 grid points in the interval [a, 6], 

so that 

Xi = a + ih, i = 0, 1, . . . , n, xq = a, x n = b, h — 



n 



We can deduce the following 

) (x)=J2C J B' J (x) 



n-l 
J=-3 

J=-3 



(4.8) 



u< 2 >(a; 



Theorem 4.3. 7/7/ze assumed approximate solution of the problem (4.5) is (4.7), then the dis- 

I n— 1 

(j=-3 



crete collocation system for the determination of the unknown coefficients {Cj}™ = * 3 is given 



by 



n— 1 n— 1 

(At)p( Xi ) ^Q^ + At^,) Y, C i B 'i( x i)+ 

j=-3 j=-3 
n-l 

[a(At)h(xi,t)-l] ^2c j B j (x i ) = -Ui-(At)n(x i ,t), 0<i<n (4.9) 

j=-3 

and boundary conditions (4.3) can be written as 

n-l 

£c^(o) = o 

J=-3 
n-l 

X)c ;l -s i (i) = e(o 

i=-3 



(4.10) 



Equation (4.9) can be written the matrix form of equations 

AC = F (4.11) 
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where 



c_ 2 



Cn-2 

\C n -i j 



Also the matrix A can be written as 

o 1 

ai(x 0 ) a 2 {x 0 ) a 3 (x 0 ) 



( =3 

h 



\ 



0 

0 
0 
0 



osi(xi) a 2 (x 1 ) a 3 (x!) 



0 
0 
0 



0 
0 
0 



\ 



0 
0 
0 



0 
0 
0 



ai(x n -i) a 2 (x n - 1 ) a 3 (x n - 1 ) 0 
0 ai(x n ) a 2 (x n ) a 3 (x n ) 

0 14 1 



where 



oti(xi) = (At)p(xi) 
a 2 (xi) = (At)p(xi) 
a 3 {xi) = (At)p(xi) 



6 



(At)q(xi) 



h 



+ [a(At)h(x h t) - 1] ■ 1, 



12 

h 2 



[a(At)h(x h t) - 1] • 4, 



6 



+ (At)q(xi 



+ [a(At)h(xi,t) - 1] • 1. 



h 2 J y '* v ' \h / 

Now we have a linear system of n + 3 equations for the n + 3 unknown coefficients, namely, 
{Cj}™Z- 3 - We can obtain the coefficient of the approximate solution by solving this linear 
system by Q-R method. 



4.3.2 Non-linear time-dependent 

Let 7](x, t, u) = h(x, t) g(u), we consider the nonlinear non-homogeneous Emden-Fowler dif- 
ferential equation 

p(x) u xx + q(x) u x + a h(x, t) g(u) + fj,(x, t) = u t (4.12) 
We can use finite difference method 



p(x) 



d 2 Ui + i , , dui 



dx 2 



+ q(x) + a h(x, t) g(u i+1 ) + fx(x, t) = Ui+1 Ui (4.13) 

ax At 
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Theorem 4.4. If the assumed approximate solution of the problem (4.13) is (4.7), then the 
discrete collocation system for the determination of the unknown coefficients is given by 



n—1 n—1 / n—1 



p(xk) ^2 Cj B"(x k )+q(x k ) ^ c j B' j (x k )+ah(x k ,ti +1 ) g ( ^ Cj Bj(x k ) )+^x k ,t i+1 ) 

J=-3 

Cj Bj(x k ) - Ui 



j=-3 j=-3 \j=-3 

yn-1 



At 

and boundary conditions (4.3) can be written as 



fc = 0,l,....n (4.14) 



n-l 



cjB'jix) = 0, for x = 0' 



'"" 3 (4.15) 

n—1 

Y CjBjix) = £(t), for x = l. 

J=-3 

Now we have a nonlinear system of n + 3 equations in the n + 3 unknown coefficients. We can 
obtain the coefficients of the approximate solution by solving this nonlinear system by Newton 's 
method. 



4.4 Examples and comparisons 

Example 4.1. [67, 64] Consider the linear problem 

2 

Ut = Uxx H — u x — (6 + Ax 2 — cos t) u, 0 < x < 1, t > 0 

Ob 

subject to boundary conditions 

« x (0,f) = 0, u(l,t) = e 1+sint , 

and the initial condition 

u(x, 0) = e x2 . 

whose exact solution is 

u(x,t) = e x2+sint . 

The maximum absolute errors at different n and different time t with At = 0.001 is given in 
Table(4.1). The numerical results for B-spline solution, exact solution and absolute error are 
illustrated in Figure 4.1 
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n 


T=0.1 


T=0.5 


T=1.0 


10 


0.002756 


0.00493 


0.007185 


20 


6.556 E-04 


0.00119 


0.001819 


30 


2.693 E-04 


5.124 E-04 


8.373 E-04 


40 


1.348 E-04 


2.742 E-04 


4.947 E-04 


50 


7.304 E-05 


1.643 E-04 


3.364 E-04 


60 


3.999 E-05 


1.049 E-04 


2.506 E-04 


70 


2.074 E-05 


6.951 E-05 


1.990 E-04 


80 


1.261 E-05 


4.683 E-04 


1.657 E-04 



Table 4.1: Results for example 4. 1 . 

Example 4.2. [67, 64] Consider the linear problem 

2 

ut + (6 - 5x 2 - Ax A ) = u xx + -u x - (5 + 4a; 2 ) u, 0 < x < 1, t > 0 
subject to boundary conditions 

u x (0,t)=0, u(l,t) = 1 + e 1+t , 

and the initial condition 

u(x,0) = x 2 + e x \ 

whose exact solution is 

u(x,t) = e x2+t + x 2 . 

The maximum absolute errors at different n and different time t with At = 0.001 is given in 
Table (4.2). The numerical results for B-spline solution, exact solution and absolute error are 
illustrated in Figure 4.2 

Example 4.3. [67, 64] Now we turn to a singular nonlinear problem 

u t = u xx + -u x - (24 i + 16t 2 x 2 ) e u -2x 2 e u/2 , 0 < x < 1, * > 0 
x v ' 

subject to boundary conditions 

u x (0,t) = 0, u(l,t) = -2 ln(l + t), 
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approximate solution 



Exact solution 




0 0 




(a) Approx. solution 



Absolute error 




0 0 



(c) Absolute error 



(b) Exact solution 



Figure 4.1: Results for example 4.1 at At = 0.01 and n = 80 



and the initial condition 



u(x,0) = 0. 



whose exact solution is 



U[X 



t) = -2 In (1 + tx 2 ) . 



The maximum absolute errors at different n and different time t with At = 0.001 is given in 
Table (4.3). 
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n 


T=0.1 


T=0.5 


T=1.0 


10 


0.002756 


0.00498 


0.00821 


20 


6.528 E-04 


0.00117 


0.00192 


30 


2.664 E-04 


4.729 E-04 


7.799 E-04 


40 


1.319 E-04 


2.309 E-04 


3.807 E-04 


50 


7.021 E-05 


1.203 E-04 


1.984 E-04 


60 


3.727 E-05 


6.203 E-05 


1.022 E-04 


70 


1.823 E-05 


2.906 E-05 


4.791 E-05 


80 


1.541 E-05 


2.431 E-05 


4.008 E-05 



Table 4.2: Results for example 4.2. 



n 


T=0.01 


T=0.1 


T=0.2 


10 


6.597 E-05 


5.805 E-04 


0.001176 


20 


1.494 E-05 


1.340 E-04 


2.807 E-04 


30 


6.481 E-06 


5.202 E-05 


1.148 E-04 



Table 4.3: Results for example 4.3. 
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approximate solution 



Exact solution 





(a) Approx. solution 



(b) Exact solution 



Absolute error 



x 10 



1.5 > 



0.5- 




(c) Absolute error 



Figure 4.2: Results for example 4.2 at At = 0.01 and n = 80 
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Chapter 5 



SOLUTION OF FOURTH ORDER TIME 



DEPENDENT PROBLEM 



This chapter is concerned with the solution based on the B-spline method of fourth-order partial 

equation with the variable coefficient of the form 

l-a\d 2 u (l + a\du . d A u fl + a\d 2 u r . , 

— )w + (— )di + ^ = (— )to* + n x >*> 0<x<1 > t>0 - 



(5.1) 



subject to the boundary conditions 



and the initial conditions 



u(0,t) = u(l,t) = 0 

du(0,t) _ du(l,t) 
dx dx 



u(x,0) = 0, 0 < x < 1, 



(5.2) 



(1 _ ff) 0u^O) = {1 _ a)go{xl 0<X< 1 



(5.3) 



where a E { — 1,1} and fx(x) > 0 is the ratio of flexural rigidity [68] of the beam to its mass per 
unit length. The case a = — 1 corresponds to the Euler-Bernoulli beam problem in elasticity. 
For more details about the formulation of this equation, see Gorman [69]. The case a = 1 
occurs in many physical models such as the theory of phase-transitions [70], shallow water 
waves [71] and an extensive literature on this subject exists [72, 73, 74, 75]. 

Many numerical methods have been developed for solving the fourth-order partial differ- 
ential equation (5.1) numerically. These methods include: finite difference methods [76], He's 
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variational iteration method [77],the Alternating Group Explicit (AGE) methods [78], Adomian 
decomposition method [79] 

The organization of this chapter is as follows. Section (5.1) contains notation, definitions 
and some results of the B-spline. The solution method is introduced in Sections (5.2) and (5.3) 
where four main theorems are presented. In Section (5.2), the error analysis of the method is 
described in detail. Section (5.3) is devoted to derivation of the discrete system. Finally, some 
numerical examples are presented in Section (5.4), where the exact solution is explicitly given. 



5.1 Description of quintic B-spline functions 



The theory of spline functions is a very active field of approximation theory and boundary value 
problems (BVPs), when numerical aspects are considered. In this paper we will be interested 
in the fifth-degree B-splines. 

Consider equally spaced knots of a partition n : a — x 0 < x\ < ■ • ■ < x n — b on [a, b]. let 
S 5 [ti] be the space of continuously differentiable, piecewise, fifth-degree polynomials on 7r.that 
is, S 5 [tt] is the space of fifth-degree splines on n Consider the B-splines basis inS^Tr]. The 
fifth-degree B-splines are defined as [50] 

x 5 , 0 < x < h 

-5x 5 + 30hx A - 60h 2 x 3 + 60/iV - 30h 4 x + 6h 5 , h < x < 2h 

10x 5 - 120hx 4 + 540/i V - 1140/i V + 1170/i 4 :r - 474/i 5 , 2h < x < 3h 

-10x 5 + \mhx A - 1260/iV 5 + 4260/i 3 x 2 - 6930/i 4 x + 438/i 5 , 3h < x < Ah 

5x 5 - 120/ia; 4 + 1140/iV - 5340/iV + 1227/i 4 x + 1097/i 5 , Ah < x < 5h 

-x 5 + 30hx 4 - 360/iV + 21Q00h 3 x 2 - QA80h 4 x + 7776/r\ 5h < x < Qh 

(5.4) 

and Bi_i(x) = B 0 [x — (i — l)h]. A detailed description of B-spline functions generated by 
subdivision can be found in[ 11]. We can simply check that each of the function B^x) G C A [a,b] 



B 0 {x) 



120/i 5 



To solve fourth-order boundary value problems B, B', B", B^\ B^ evaluated at the nodal 
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points are needed 



Bi( Xj ) = < 



66, 
26, 

1, 
0, 



similarly, we can show that: 



o, 


if; 


i = j 




_i_ 50 

^ h ' 


if; 


i-j 


= ±1 




if; 


i-j 


= ±2 


o, 


if; 


i-j 


= ±3 


o, 


if; 


i 


— 3 



if; 

if; 
if; 
if; 



_i_ 120 

31 fe3 , 



_i_ 60 

0, 



if; i-j = ±l 
if; i- J = ±2 

if; i-j = ±3 



i = 3 

i-3 = ±1 
i-j = ±2 

i-3 = ±3 



120 

fe2 > 



40 

h 2 ' 

20 
ft 2 ' 



5} 4) (*i) = < 



(5.5) 



720 

h 4 ' 



«/; i = j 

if; i-j = ±l 
if; i- j = ±2 

if; i-j = ±3 
(5.6) 

if; i = j 



120 

ft 4 ' 



if; i-j = ±2 

if; i-j = ±3 
(5.7) 



5.2 Error analysis 
5.2.1 Governing equation 



Setting 



d 2 u u i+ i - 2ui + Ui-i 



and 



du 
~dt 



dt 2 (At) 2 
Equation (5.1) may be approximated by the following scheme 



Uj - Uj-i 

At 
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'1-a 



(At) 2 



+ 



l + <7 



At 



d A u i+l 



1 + a\ d 2 u 
2 J ~I^ 2 



(5.8) 



i+l 



d A u 



1 + er\ gPm 



' + ( TTT^TTTo ) = l/>(x,t i+1 ) 



where 

= f(x,t i+1 ) 



1 - (7 

2(Atp 



l + o- 
2At 



Mi - 



A2(At) 2 



l + cr 
2At 



(5.9) 



Ui-l 



Suppose S(x) is the quintic B-spline interpolating u(xk, t i+ i) at the % + 1 time level then 



n-1 



S'(x) = ^ c kB 5 ,k( x ) 



(5.10) 



fe=-5 



5.2.2 Truncation error of B-spline function 

Lemma 5.1. Let's denote by Ui = u(xi) and S^ = S&\xi) for all p where S p = D P S and 
u(xi) S'(xi). Further, we define A by AS* = Sj_ 2 + 26Sj_i + 66Sj + 26S i+ i + Sj+2 for any 
function S evaluated at the nodes n then we have the following recursive relations. 



and 



AS" = i [5w(xj_ 2 ) + 50w(a; i _i) - 50u(z i+ i) - u(x i+2 )\ , 
20 

AS"' = — [u(xi- 2 ) + 2w(xi_i) - 6u(xi) + 2u(x i+ i) + w(x i+2 )] , 
/// 60 

AS"" = ^ [w(xj_ 2 ) - 2u(xi_i) + 2u(x i+1 ) - u(x i+2 )\ , 



120 

AlS<(4) = [ u ( x »-2) _ 4u(xi_i) + 6u(xj) - 4u(x i+ i) + u(x i+2 )] 



(5.11) 
(5.12) 
(5.13) 

(5.14) 



Theorem 5.1. Le? S be the unique quintic spline for given function u , then for uniform parti- 
tions, we have for i = 0(l)n then 



S$ = uj + 0 (Zi 8 ) , 



s'! = u'j + 



720 



-u 



(6) 



+ o (h 6 ) , 



qlll _ III „,(7) 



//, ' +0(/i 6 ) , 



s. 



( 4 ) _ „.(4) 



/l 2 



,(6) 



(5.15) 
(5.16) 
(5.17) 
(5.18) 
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Proof. It is enough to prove the equation (5.15) since the proofs of the remaining cases are 
similar. By using the operator notation E(S(xi)) = S(x i+ i) and the above relation (5.11) may 
be written as [12] 



[E- 2 + 26 E- 1 + 66 + 26 E + E 2 ] S'fa) = (J^j [E~ 2 + 10 E' 1 - 10 E - E 2 ] 



U{Xi) 

(5.19) 



since E = e hD where D = -f- substitute in the above relation 

ax 



[e~ 2hD + 26 e~ hD + 66 + 26 e hD + e 2hD ] S'( Xi ) = [e~ 2hD + lOe"^ - 10e hD - e 2hD ] u(x t ) 

(5.20) 

Expand in Maclaurine series in power of hD then adding similar terms we will get 



S'(xi) 



2Ah 



AO 

2AhD + 6h 3 D 3 + ^-h 5 D 5 + —h 7 D 7 



42. 



7! 



l+(I^+^ + 



2Ah 



42 48 
2AhD + Qh 3 D 3 + —h 5 D 5 + —h 7 D 7 
60 7! 



1-1^ + ^ + 



+ ( h 2 D 2 + —h 4 D 4 + 
\4 720 



+ 



u(xi) 



where Du(x) = u x (x) and _D 7 -u = u X xxxxxx(£) and so on 

5237 



s , w = ( D +S k ' DT+ -)" (l,) 



which give directly (5.15). 



(5.21) 

□ 



Theorem 5.2. 5y wsmg the combination of the finite difference approximation and the quintic 
B-spline method, the truncation error for the time-dependent equation (5.1) is 

O (h 2 + (At)) when a = 1 
O (h 2 + (At) 2 ) when a = -1 
Proof. By applying the finite difference scheme with time step At to equation (5.1) 

A [u(U + At, x) - 2u(U, x) + u(U - At, x)} + ( j [u(U, x) -u(U- At, x)} 



2(Aty 



2At 



+M*)« (4) = ( ^) « (2) + /OMi) 



(5.22) 
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when apply Taylor expansion in the above relation 



l-<7 

2(Atj 2 



(At) 2 u tt + y^umt + ■■■ + 0(At 6 ) 



12 



+ 



1 + 0" 



(Aty 

3! 



u ttt ■■■ + 0(At 4 



2At 

+/i(x)M (4) 



(At) 2 

* ~ 2T Utt+ 



collecting similar terms together 



1 - a \ ( 1 + a \ , , M (I — a 

U tt + [ - ) «t + ; + 



(At) s 
12 



««« + • • • + 0(At 4 ) 



+ 



l + [ At (At) 2 . A „, 



1 + 0- 



u w + /(x,ti) 
(5.24) 



then use the relations (5.11-5.14) and substitute in (5.24) 



l-a\ fl + cr\ (\-a 

u tt + —z— u t + 



2 

1 + 0 
2 

H(x) 



12 



utttt + ■■■ + 0(At 4 ) 



+ 



At (At) 2 . A „. 



+ 



1 + 0 



u" + A_^( 6 ) + 0(/> 6 ) 



(5.25) 



and the error can be written as 



e = 



1 + 0 



+ 0(/i e ) 



720 



1-0 



-/i(x) 

(At) 2 

12 
1 + a 



+ • • • + 0(At 4 ) 



At (At) 2 n , A 

2r«« + ^p M *« + • • • + °(At 3 ) 



(5.26) 



when we substitute about a = 1 this give order of error O (h 2 , (At)) and when a = — 1 this 



give order of error O (h 2 , (At) 2 ) and this complete the proof. 



□ 



5.3 Scheme of the method 



Let the solution u(x, t i+ i) of the problem (5.9) be approximated by 



u 



n-l 
j=-5 



(5.27) 
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where Cj are unknown real coefficients and Bj(x) are the fifth-degree B-spline functions. Let 
xoi xi, x n be n + 1 grid points in the interval [0, 1] so that x k = k h, k = 0, 1, 2, n where 

x 0 = 0, x n = 1 and h — — 

n 

It is required that the approximate solution satisfies the differential equation at the points x = 
x k . and we can easily deduce the following 

n— 1 n— 1 

"'(Zfc) = 5^ C j B j( X k), u"(x k ) = E CjB"(x k ), 

j= ~ 5 S= ~*, (5.28) 

n— 1 n— 1 

« (3) (^) = E c ^f « (4) w = E c 3 B f\^). 

j=-5 j=-5 

Theorem 5.3. If the assumed approximate solution of the problem (5.9) is (5.27)then the dis- 
crete collocation system for the determination of the unknown coefficients {Cj}™~^ 5 is given 
by 

n—l / ^ \ n— 1 /• j \ n— 1 

j=—5 • j=— 5 J=— 5 

(5.29) 

and boundary conditions (5.2) can be written as 

n— 1 n— 1 

E CjBj(x) = 0, E c i 5 j(» = °> for i = 0 

(5.30) 

n— 1 n—l 

E Cj-Sj(a;) = 0, E 9 5 j( ;r ) = °' for x = 1 

j=-5 j=-5 

Proof. We replace each term of (5.9) with its corresponding approximation given by (5.28) and 
substituting x = x k and applying the collocation to it. □ 

From the initial condition from i = 0, we have 

u Q = u(x, 0) 

Mi = (At) /jL 2 (x) + Uq 

Then the system in (5.30) takes the matrix form 

AC = Q (5.31) 
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where 



C = 



/ \ 

C-5 
C_ 4 



Cn-2 
C n -1 



Q = 



o 
o 



also the matrix A can be written as 



/ 



A = 



26 

_50 
h 



66 
0 



26 

50 
h 



5 
h 



ip(x 0 ,t i+1 ) 

tp(x n ,t i+1 )) 
0 
0 

0 0 
0 0 



ai(x 0 ) «2(^o) «3(^o) «4(^o) «5(^o) 0 



0 



0 
0 
0 



0 
0 
0 



0 
0 
0 



0 
0 
0 



0 
0 
0 



0 
0 
0 



0 
0 
0 



0 ai(x n ) a 2 (x n ) a 3 (x n ) a A (x n ) a 5 (x n ) 



0 
0 



50 
' h 



26 



0 

66 



50 
h 

26 



/ 



Also the coefficients of the matrix A have the following form 



ai(x k ) 


= K x k) 




= n{x k ) 


a 3 (x k ) 


= n{x k ) 


a 4 (x k ) 


= n(x k ) 


0-s(x k ) 


= K x k) 



120\ 

/-480 

/720\ 

/-480 
V h 4 

m 



l + a 
2 
1 



a 



2 

l + a 
2 

l + a 
2 

l + a 



40 \ 
-120' 



h 2 
40 



+ 



1-a 
2(Atj2 

(— 
\2(At) 2 

1 - a 



2(At) 2 
1 -a 



— 

h 2 ) + \ 2(At) 2 



.(26), 
) -(66), 
.(26), 

(!)■ 



(5.32) 



Theorem 5.4. The matrix A in the system (5.31 )is nonsingular. 



Proof. The basis function B,i(x) is defined only in the interval [xj_ 2 , x i+2 \ and outside of 
this interval it is zero. Therefore,!?; (x) is having non- vanishing values at the mesh points 
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{xi-2,Xi-i, %i, Xi + i,Xi +2 } and at other mesh points the value of Bi(x) is zero. It is clear that 
from the definition of quintic B-splines defined in section 2, the value of Bi(x) at Xj is domi- 
nating for j = i when compared with the values of B^x) when j ^ i. The derivatives of B^x) 
up to fourth order also have the same nature at the mesh points as in the case of Bj(x) Using 
these facts, we can say that the matrix A is a tri-diagonal band matrix with nonzero entries and 
dominant principal diagonal elements. Hence the matrix is nonsingular. □ 

Now we have a linear system of n + 5 equations for the n + 5 unknown coefficients, namely, 
{Cj}"~^ 5 . We can obtain the coefficients of the approximate solution by solving this linear 
system by Q-R method. 

5.4 Examples and comparisons 

In this section we introduce 3 examples and for each example the error is tabulated and figures 
of Approximate, Exact and Error are introduced. 

Example 5.1. In the case, a = 1 and //(#) = e x equation (5.1) becomes 



u t + e u. 



xxxx 



u xx + f(x,t) 



where 



f(x,t)= 3t 2 (l - t) (. 



x — X 




subject to the boundary conditions 



u(0,t) = 0 



u(l,t) = 0 



u x (0,t) = 0 



and the initial condition 



u(x,0) = 0 



whose exact solution is 
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The maximum absolute error is displayed in Table (5.1) for different n at time 1.0 with time 
step At = 0.001. The numerical results for B-spline solution, exact solution and absolute error 
are illustrated in Figure (5.1) 



n 


Maximum Absolute Error 


10 


7.66873 E-05 


20 


1.91713 E-05 


30 


8.52628 E-06 


50 


3.07671 E-06 


80 


1.20885 E-06 


100 


7.77820 E-07 



Table 5.1: Results for example 5.1. 



Example 5.2. [80] In the case, a = — 1 and ji(x) — 1 + x equation (5.1) becomes 

u tt + (l + x)u xxxx = f(x,t) 

where 

f(x,t) = -72(1 + a;) (l - 5x + 5x 2 ) t sint + [x - x 2 ) 3 (2 cos t - t sin t) 
subject to the boundary conditions 

u(0,t) = u(l,t) = 0 
u x (0,t) =u x (l,t) = 0 

and the initial condition 

u(x, 0) = u t (x, 0) = 0 

whose exact solution is 

u(x, t) — (x — x 2 ) 3 t sin t 

The maximum absolute error is displayed in Table (5.2) at time 1.0 with time step At = 0.001. 
The numerical results for B-spline solution, exact solution and absolute error are illustrated in 
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approximate solution for Ex(1) 



Exact solution for Ex(1) 




0 0 




(a) Approx. solution 



(b) Exact solution 



Absolute error for Ex(1) 



1.5- 




0.5- 



(c) Absolute error 



Figure 5.1: Results for example 5.1 at At = 0.01 and n = 100 



Figure (5.2) 



Example 5.3. [81] in the case, a = —I and jj,{x) = 1 + x equation (5.1) becomes 



u tt + (1 + x)u xxxx = f{x,t) 



if 



f(.r. i ) = [ x 4 + x 3 — —x 7 ) cos / 
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n 


Maximum Absolute Error 


10 


0.0013128 


20 


3.2854E-04 


40 


8.2418E-05 


70 


2.7162E-05 


120 


9.4873E-06 



Table 5.2: Results for example 5.2 
subject to boundary conditions 

u(0,t) = 0 

«x(0,t) = 0 Ux (l,t) 



, 6 7 
u(l, t) = —x cost 

cost 



120 ' 



and the initial condition 



u(x,0) = -x 7 



then the exact solution is 



u t (x,0) = 0 



1 



u(x,t) = ——x 7 cost 
v ' ; 840 



the maximum absolute error at time 1.0 with time step At = 0.001 are displayed in table (5.3) at 
time 1.0 with time step At = 0.001. The numerical results for B-spline solution, exact solution 
and absolute error are illustrated in Figure. 5.3 



n 


Absolute error 


10 


6.22286 E -006 


20 


1.46198 E-006 


30 


6.42049 E-007 


80 


9.01539 E-008 


100 


5.79109 E-008 



Table 5.3: Results for example 5.3. 
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approximate solution for Ex(2) Exact solution for Ex(2) 




(a) Approx. solution (b) Exact solution 

Absolute error for Ex(2) 



-5 

x 10 




(c) Absolute error 

Figure 5.2: Results for example 5.2 at At = 0.01 and n = 100 
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(c) Absolute error 



Figure 5.3: Results for example 5.3 at At = 0.01 and n = 100 
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Chapter 6 



NUMERICAL SOLUTIONS FOR 



FALKNER-SKAN EQUATION 



6.1 Introduction 

Finding the numerical solution of nonlinear BVPs on infinite intervals is one of the problems at- 
tracting many scientists. The nonlinear third-order Falkner-Skan equation is a famous example 
of the BVPs on infinite intervals arisen in many branches of sciences, e.g. applied mathematics, 
physics, fluid dynamics and biology. We are to solve the following Falkner-Skan equation 



1 - 



dm 

dr) 



0, ?]G[0,oo) (6.1) 



drf drf 
subject to boundary conditions 

/(0) = /'(0) = 0 (6.2) 

f'(rj) = 1 as i] ->■ oo (6.3) 

Which arises in the study of viscous flow past a wedge of angle f3 n and f3 0 and f3 are constants. 

In the case of ft — 0, the Blasius equation is obtained, this equation is perhaps one of the 
most famous equations of fluid dynamics and represents the problem of an incompressible fluid 
that passes on a semi-infinity flat plate. In the case of accelerating flows f3 > 0 , the velocity 
profiles have no points of inflection, whereas in the case of decelerated flows (3 < 0. Physically 
relevant solutions exist only for —0.19884 < j3 < 2. 



70 



CHAPTER 6. NUMERICAL SOLUTIONS FOR FALKNER-SKAN EQUATION 



Approximate solutions of the nonlinear third order boundary value problem (6.1)-(6.3) are 
obtained by many researchers, see for example,[82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92] and 
the references there in. 

The organization of the chapter is as follows: Section (6.2) contains notation, definitions 
and some results of the quartic B-spline. In Section (6.3), we summarize the application of B- 
spline method for solving Falkner-Skan equation and in Section (6.4) compare it with existing 
methods in the literature like [83, 84, 93, 94]. 



6.2 Description of quartic B-spline functions 

We first use basic assumptions to derive the expressions for the quartic uniform B-splines di- 
rectly and without mentioning knots. We then show how to extend the derivations to uniform 
B-splines of any order. Following this, we discuss a different, recursive formulation of the 
weight functions of the uniform B-splines. 

If the interval [a, b] is partitioned into n + 1 uniformly spaced points x rn such that a = x 0 < 

xi < x 2 ... < x n = b and h = The quartic B-splines B m (x),m = —2, —3, , n + 1 at the 

knots x m are defined as: 



B m (x) 



h 4 



x - x m - 2 ) A if: 

x - x m - 2 f - 5(x - x m _i) 4 , if: 

x - x m ^ 2 ) A - 5(x - rc m _i) 4 + 10(x - Xm) 4 if: 

x m+3 - x) 4 - 5(x m+2 - x) 4 if: 

x m+3 - x) 4 if: 

0 if: 



[^m— 2i X m — l] 

[x m , X m +i\ 

[x m +l, X m+2 ] 
[Xm+2-, Xm+3\ 

otherwise 



(6.4) 



and the set {B^ 2 , B_i, B n+1 } of quartic B-splines forms a basis over the interval [a, b]. 
The values of quartic B-spline function and its derivatives are summerized in table (6.1) 
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X 


%m-2 


•Em—I 






x m+2 




B 


0 


1 


11 


11 


1 


0 


B' 


0 


4 
h 


12 
h 


12 
h 


4 
h 


0 


B" 


0 


12 
h 2 


12 
h 2 


12 

h 2 


12 

h 2 


0 


B'" 


0 


24 
h 3 


72 

h 3 


72 
h 3 


24 
ft 3 


0 



Table 6.1: The values of quartic B-spline and its derivatives at nodal points . 

6.3 Description of B-spline functions 
6.3.1 The derivation of Falkner-Skan equation 



In this section, we introduce two-dimensional steady- state boundary layer equation with 
boundary conditions as follows 

In two-dimensional, the equation of continuity and the laminar boundary layer equations for 
the steady flow of an incompressible viscous fluid over a wedge are given by [95] 

dv dw 
ox oy 

dv dv dll d 2 v 
ox oy ax oy z 

where v; w are components of velocity in x and y direction of the fluid flow, /i is the viscosity, 
U (x) is the velocity at the edge of the boundary layer. We consider a general case of a power 
law free stream velocity, that is, U (x) = [Too (f )"\ where [Too is uniform free stream velocity, 
L is the length of the wedge, x is measured from the tip of the wedge and m is the Falkner-Skan 
power-law parameter. The boundary conditions are given by 

v(x, 0) = w(x, 0) = 0 

(6.7) 

v(x,y) — ?■ U{x) as y — > oo 

The continuity equation (6.5) is automatically satisfied by the stream function u(x, y) such that 

dv dv 
v = — , and w = -— (6.8) 

dy dx 
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and the momentum equation (6.6) takes the form 

dv d 2 v dv d 2 v 



dU d 3 u 
dy dy dx dx dy 2 dx ^ dy 3 



(6.9) 



By the similarity transformation 

f(v) = 



1 +771 L r 



fj,U 0 



1+m 
X 2 



7} 



1+m U 0 



y 



1 — m 
X 2 



(6.10) 



2 n L r ' 

where / is a dimensionless stream function and rj is a dimensionless distance from the edge, 
called similarity variable, the partial differential equation (6.9) reduces to third order nonlinear 
ordinary differential equation (6.1) where (5 n is the wedge angle and is related to the Falkner- 
Skan power-law parameter m through the expression (5 = and prime ' denotes differentia- 
tion with respect to 77. 

The equation (6.1) is well known Falkner-Skan boundary layer equation [95] in which f'(ri) 
defines the dimensionless velocity component indirection and f"(rf) defines the dimensionless 
shear stress in the boundary layer. The boundary conditions (6.7) can be written as (6.2)-(6.3). 



6.3.2 Transformation to a finite domain: 



First of all, we replace the boundary condition at infinite (6.3) with a free boundary condition 
as Asaithambi did in [82]: 

df 1 

— = 1 at 77 = 7700 

d 2 f (6 - U) 

= 0, at 77 = 77^ 

dr] z 

in which 77^ is the unknown free boundary (truncated boundary). Then the original problem 
(6.1)-(6.3) becomes the free boundary problem (6. l)-(6.2) and (6.11)defined on a finite interval, 
where 77^ is to be determined as a part of the solution. We add another initial value condition 
at 77 = 0 

d 2 f 

-t4 = «. at ^ = ° ( 6 - 12 ) 

df] z 

The coordinate transformation adopted here is as follows: 

c = ^ 
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This maps the physical domain [0, r^] to the fixed computational domain [0, 1]. If we let ip = 
then after some algebraic manipulation, equation (6.1) can be rewritten: 



^ + ^LM% + Vim - (p 2 ] = 0 (6.13) 
the boundary conditions are then: 

^(0=^(0 = 0 at C = 0 (6.14) 

and conditions (6.1 1) will be 



^'(0 = 1 at C = l, 

^'(C) = 0 at C = l- 



The additional condition will be 

6.3.3 The numerical method 



r^oo a, at ( = 0 



(6.15) 



Let 

n+l 

lKC0 = ^CjBjiCi) (6-16) 

be an approximate solution of equation (6.13), where Cj are unknown real coefficients and 
Bj(() are fourth-degree B-spline functions. 

Let Co, Ci) • • • , Cn be n + 1 grid points in the interval [0, 1], so that 

Q = ih, i = 0,l,...,n, Co = 0, ( n = 1, /i 



And we can easily deduce 



n+l 



J=-2 
n+l 



J=-2 
n+l 

J=-2 
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The nodal values B, B', B" and B'" at the knots Q are derived from table (6.1) Substitute from 
(6.16-6.17) into (6.13) yields. 



n+1 / n+1 \ / n+1 

E cjBj'itiWoofo E c > B i&) E ^-(co 

j=-2 \j=-2 J \j=-2 



n+1 

1 - ( E c > B 'M> 

J=-2 



(6.18) 



And boundary conditions will be 

n+1 



Ec^(0) = o, 

J=-2 
n+1 

E Cj - o- 

J=-2 
n+1 

EQ^(0) = o, 

J=-2 



i = 0, 
% = 0, 
i = n. 



(6.19) 



Now we have a nonlinear system of n + 5 equations in the ra + 4 unknown coefficients and the 
unknown value r]^. We can obtain the coefficients of the approximate solution by solving this 
nonlinear system by the following iterative process. 

6.3.4 Solution of Nonlinear system 

Assuming that ^ is available then (6. 18)-(6. 19) will be used to determine the (iV+4) unknown 
coefficient { Cj } We begin by letting. [80, 83] 

Fl(c_ 2 ,C_i, ...,C n+ i) 0 



F(C) = 



F 2 (c_ 2 ,c_i, ...,C n+ i) 



Fn+i{C-2, C-l, C n+ i) 



V 



/ 



(6.20) 



V / 



Where C is the column vector of independent variables and F is the column vector of the func- 
tions Fj, with Fj(C) = Fj(c_ 2 , c_i, c n+1 ), 1 < i < n + 4. The number of functions that are 
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set to zeros is equal to the number of independent variables. A very good method for solving 
equation (6.20) is Newton method. 

Let c^' fc ) be the guess at the solution for iteration j and F^' fc ^ denote the value of F at the j th 
iteration. Assuming that || F^ ,fc ^ || is not too small, we seek update vectors Ac^'^ 



c -2 



c (j+l,fe) _ c (j,k) ^ c (j,k) 



such that F(c J ' +1 ' fc ) = 0. Using the multidimensional extension of Taylor's theorem to approxi- 
mate the variation of F(c) in the neighborhood of c l gives 

F(c j ' k + Ac 1 *) = F(c (j) ' fe ) + F'(c (j ' ) > fc )Ac (j ' ) ' fc + o(\\ Ac (j) ' k f) (6.22) 



c -2 



J* 



e> 



,k 



f Ac^ 2 ^ 



+ 



Ad:' 



(6.21) 



where F'(c^' fe ) is the jacobian of the system of equations 



/ dFi 



F'(c)=J(c) 



<9c_ 



dF 2 

dc-o 



W £( c ) 



OF, 



dF± 

dc. 



n+l 



<9F 2 



9c 



n+l 



(c) 



(c) 



\ 



<9c n 



+i 



(6.23) 



V / 
Neglecting higher order terms and designating as the Jacobian evaluated at c^' k \ We 

can rearrange equation(6.22). 



F( c ^> fe + Ac U) ' k ) = F(c {j) ' k ) + J ij) Ac (i 



(j),k 



(6.24) 
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The goal of Newton iterations is to make F(c^' fe + Ac^' k ) = 0 , so setting that term to zero in 
the preceding equation gives 

j(i) Ac 0'),fc = _F( c W' fe ) (6.25) 

Equation (6.25) is system of Q linear equations in the Q unknown Ac^' k . Each Newton iter- 
ation step involves evaluation of the vector the matrix and the solution to equation 
(6.25). A common numerical practice is to stop the Newton iteration whenever the distance 
between two iterates is less than a given tolerance; i.e., when || c (j+1 )' fc — c^ ,fe ||< e\. 
Let [83] 

v ( VoD (k)) = ^"(l) = 0 (6.26) 
by applying finite difference it will converted to. 

V ( Voo (k)) = <+i 2 ;; + ^ (6.27) 



then apply successive iterations to determine rj^ using secant method 

v£ +1) = v<$ - v^l fc) ) 



(k) (fe-1) 



(6.28) 



for suitable initial guesses 77^ and rj^ . A common numerical practice is to stop the secant 
iteration whenever the distance between two iterates is less than a given tolerance; i.e., when 

II (<+l) (0 IK 

II '/oo '/oo || ^ C2- 



6.3.5 Algorithm 



We are now ready to present the algorithm of our method: 



initialize rj^ = r]^\ 



/oo , 



F(c j ,r] 00 ) = ¥(cj,r] ( 0 
initialize c = c^°\ 
fori = 0, 1,2, •••(p(rj 00 ) = 0, 
for A; = 0,1,2, ■■■¥ {k) <\c j ) = 0, 
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• if || F (fc),i || is small enough, stop, 

• compute J (fc) '* 

• solve jW-'AcW-* = -F(c^), 

• c (fe+l),i _ c (k),i _|_ ^ c (k),i 

• compute ipi 

(i) (i-l) 
vfoS)-^- 15 ). 

• end 

• end 



6.4 Results and discussion 



In this section we will introduce the solution Falkner-Skan problem for multiple values of f3 
according to the values stated in section 6.1 and for f3 0 — 1. Where the solution is physically 
relevant —0.19884 < (5 < 2. We compare our results with other methods stated above. And the 
value of a is calculated and tabulated in table (6.4) where a — ^4 but in order to calculate a in 
terms of ip the it will be approximated using finite difference as following. 



a 



j) 0 - 2^1 + ^ 2 
h 2 



(6.29) 
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Figure 6.1: The velocity profile corresponding to (3 G {0, 0.5, 1, 2} 




Figure 6.2: The velocity profile corresponding to (3 e {—0.1, —0.12, —0.15, —0.18} 
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p 


B -spline 


method 
[83] 


method 
[84] 


Chebyshev 
[93] 


ADM [94] 


2.0000 


1.687210 


1.687222 


1.687222 


1.687218 


1.683050 


1.0000 


1.232587 


1.232589 


1.232588 


1.232588 


1.196315 


0.5000 


0.9276795 


0.927682 


0.927680 


0.927680 


0.807826 


0.0000 


0.469600 


0.469601 


0.469600 


0.469601 


0.491960 


-0.1000 


0.319270 


0.319270 


0.319270 


0.319272 


0.320189 


-0.1200 


0.281761 


0.281762 


0.281761 


0.281765 


0.278509 


-0.1500 


0.216361 


0.216360 


0.216361 


0.216376 


0.224771 


-0.1800 


0.128636 


0.128637 


0.128636 


0.128682 


0.129001 


-0.1988 


0.005237 


0.005217 


0.005220 


0.005270 


0.005227 



Table 6.2: Comparison of results for a of Falkner-Skan problem. 
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CONCLUSION AND FUTURE WORK 

7.1 Conclusion 

In this study, we compared the performance of the collocation method using B-spline bases and 
other methods for solving boundary value problems. Numerical experiments are presented. 

The B-spline collocation method is a simple method with high accuracy for solving linear 
and non-linear problems. So it may be easily applied by researchers and engineers familiar with 
B-spline function. 

B-spline method has been successfully employed to obtain the approximate solutions of 
various linear and nonlinear second and fourth order time-dependent problems and some nu- 
merical examples are presented and we introduce complete error analysis to the model of time- 
dependent equation which is stated in theorem 3.2. Moreover, the B-spline method not only 
used in a direct way but also requires less computational work in comparison with the other 
methods. 

The B-spline technique is applied to solve Emden-Fowler problem which contains singular 
term .the analysis of the method exhibits the reliable applicability of B-spline method to solve 
linear and nonlinear problems with singular feature. The difficulty in this type of equations, due 
to the existence of singular point at x = 0, is overcome here. Moreover, We have presented a 
few examples to illustrate the method and compared the solution with the exact ones to illustrate 
the effectiveness of the proposed method. We apply the quartic B-splines method to Falkner- 



82 



CHAPTER 7. CONCLUSION AND FUTURE WORK 



Skan equation and write an algorithm depends on newton and secant methods to solve the 
non-linear system of equations results after applying the B-spline collecation method to the 
equation. Moreover, We have presented a comparison of B-spline results with other methods to 
verify the accuracy of the presented technique. 

7.2 Future Work 

There are still existing many open questions related to the topics of this thesis. For instances, 

• study of integral and integro differential equations. 

• Study of partial differential equation with nonlocal boundary conditions. 

• Study of B-spline method to solve Stefan problem. 

• Study of B-spline-collocation method to solve fractional differential equation. 

• study of Sturm-Liouville eignvalue problems. 

• Study of B-spline-collocation method to solve fractional integral equation. 

• Study of B-spline-collocation method to solve the BagleyTorvik equation 

• using higher order B -splines to solve higher order ordinary differential equations 

• use non-polynomial splines to solve different types of partial differential equations 
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